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Nomenclature 



A 


Thermal flux production vector, (1.7e) 





Buoyancy parameter, g(3 


B 


Pressure-temperature-gradient-vector, (1.7f) 




Thermal expansion coefficient 


C 


Energy conversion tensor, (1.7b) 


7m 


Relaxation frequency of ry , i = j 


v/vt 


Substantial derivative, d/dt + U- V 


7m 


Relaxation frequency of , i =/= j 


D/Dt 


Mean substantial derivative, d/dt + U- V 


7rd 


Relaxation frequency of F 


E K 


Turbulent kinetic energy per unit mass, \u\ 2 /2 


luu 


Relaxation frequency of E K 


E e 


"Temperature energy" per unit mass, 2 /2 


lee 


Relaxation frequency of E e 


F 


Turbulent thermal flux per unit mass, (u8) 




Dissipation tensor of ry, (1.6) 




Thermal flux at zero elevation z = 


e 


Dissipation vector of F, (1.6) 


g 


Gravity acceleration, g = —gz 


£ 


Dissipation of E e , (1.6) 


L 


Monin-Obukhov length, / '/3F* 




Deviation of potential temperature from BRS 


i 


Outer scale of turbulence, external parameter 


e 


Mean potential temperature, (Qd) 




Rate of Reynolds stress production, (1.7a) 





Fluctuating potential temperature, 0^ — (Qd) 


P, P, P* 


Total, fluctuating and zero level pressures 


0, t 


Potential temperature 


Pr T 


Turbulent Prandtl number, v t /Xt 




at zero elevation, F*/u* 


Rifiux 


Flux Richardson number, /3F z /t xz S u 


A* 


Viscous lengthscale, v/u* 


Rigrad 


Gradient Richardson number, (3S e / S 2 


V 


Kinematic viscosity 




Pressure-rate-of-strain-tensor , ( 1 . 7c) 




Turbulent viscosity 


s v 


Mean velocity gradient, dll/dz 


p 


density of the fluid 


S e 


Mean potential temperature gradient, dO/dz 


Tij 


Reynolds stress tensor, (uiUj) 


T 


Molecular temperature 


T* 


Mechanical momentum flux 


U 


Velocity field 




at zero elevation (at the ground) 


u 


Mean velocity, (U) 


X 


Kinematic thermal conductivity 


u 


Fluctuating velocity, U — U 


Xt 


Turbulent thermal conductivity 




(Wall) friction velocity, y 7 ?. 


BRS 


Basic Reference State 


X 


horizontal (streamwise) unit vector 


z 


vertical (wall-normal) unit vector 
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Introduction 

The lower levels of the atmosphere are usually strongly 
influenced by the Earth's surface. Known as the atmo- 
spheric boundary layer, this is the part of the atmosphere 
where the surface influences the temperature, moisture, 
and velocity of the air above through the turbulent trans- 
fer of mass. 

The stability of the atmospheric boundary layer de- 
pends on the profiles of the density and the temperature 
as a function of the height above the ground. During 
normal summer days the land mass warms up and the 
temperature is higher at lower elevations. If it were not 
for the decrease in density of the air as a function of 
the height, such a situation of heating from below would 
have been always highly unstable. In fact, the boundary 
layer is considered stable as long as the temperature de- 
creases at the dry adiabatic lapse rate (T" ~ — 9.8°C per 
kilometer) throughout most of the boundary layer. With 
such a rate of cooling one balances out the decrease in 
density. With a higher degree of cooling one refers to 
the atmospheric boundary layer as unstably stratified, 
whereas with a lower degree of cooling the situation is 
stably stratified. Stably stratified boundary layer oc- 
curs typically during clear, calm nights. In extreme cases 
turbulence tends to cease, and radiational cooling from 
the surface results in a temperature that increases with 
height above the surface. 

The tendency of the atmosphere to be turbulent does 
not depend only on the rate of cooling but also on the 
mean shear in the vertical direction. The commonly used 
parameter to describe the tendency of the atmosphere 
to be turbulent is the "gradient" Richardson number 
(Richardson, 1920), defined as 

_ 0dB{z)/dz 
Rlgrad= [dU x /dz]* ' (0 } 

where x is the stream-wise direction, z is the height 
above the ground, @(z) is the mean potential temper- 
ature profile, (which differs from the mean temperature 
profile T(z) by accounting for the adiabatic cooling of the 
air during its expansion: d<d(z)/dz = dT{z)jdz + \T'\), 

(3 = (3 g is the buoyancy parameter in which f3 is the 
adiabatic thermal expansion coefficient (for an ideal gas 
(3 = 1/T), and g is the gravitational acceleration. The 
mean shear dU x /dz is defined in terms of the mean ve- 
locity U, which in the simplest case of flat geometry de- 
pends only on the vertical coordinate z. The parameter 
Rigrad represents the ratio of the generation or suppres- 
sion of turbulence by buoyant production of energy to 
the mechanical generation of energy by wind shear. 

This paper is an extended presentation of an invited 
talk given on International Conference "Turbulent Mix- 
ing and Beyond" devoted, in particular, to the problems 
of fluid dynamics, turbulence, geophysics and statistics, 
that are long-standing challenging tasks. Here we con- 
sider the description of stably stratified turbulent bound- 



ary layers (TBL), taking as an example the case of sta- 
ble thermal stratification. Since the 50's of twentieth 
century, traditional models of stratified TBL generalize 
models of unstratificd TBL, based on the budget equa- 
tions for the kinetic energy and mechanical momentum; 
see reviews of Umlauf and Burchard (2005), Weng and 
Taylor (2003). The main difficulty is that the budget 
equations are not closed; they involve turbulent fluxes of 
mechanical moments (known as the "Reynolds stress" 
tensor) and a thermal flux F (for the case of thermal 
stratification): 

7fc = <Uji*j) , F={u6), (0.2) 

where u and 9 stand for the turbulent fluctuating velocity 
and the potential temperature with zero mean. The na- 
ture of the averaging procedure behind the symbol (• • • ) 
will be specified below. 

Earlier estimates of the fluxes (0.2) are based on the 
concept of the down-gradient turbulent transport, in 
which, similarly to the case of molecular transport, a 
flux is taken proportional to the gradient of transported 
property times a corresponding (turbulent) transport co- 
efficient: 

t xz = -is T dU x /dz, u T rj C v 4\/t^ , (0.3a) 
Fx = -XTdB/dz, XtwC x 4V^I, etc. (0.3b) 

Here the turbulent-eddy viscosity v T and turbulent ther- 
mal conductivity Xt are estimated by dimensional rea- 
soning via the vertical turbulent velocity y/T^l and a scale 
l z (which in the simplest case is determined by the ele- 
vation z). The dimensionless coefficients C v and C x are 
assumed to be of the order of unity. 

This approach meets serious difficulties (Zeman, 1981), 
in particular, it predicts full suppression of turbulence 
when the stratification exceeds a critical level, for which 
Rigrad ~ 0.25. On the other hand, in observations of the 
atmospheric turbulent boundary layer turbulence exists 
for much larger values than Ri gra d = 0.25: experimen- 
tally above Ri gra d = 10 and even more (see Galperin 
et al. (2007) and references therein). In models for 
weather predictions this problem is "fixed" by introduc- 
ing fit functions C„(Ri gra d) and C x (Ri gra d) instead of the 
constant C v and C x in the model parametrization (0.3). 
This technical "solution" is not based on any physical 
derivation and just masks the shortcomings of the model. 
To really solve the problem one has to understand its 
physical origin, even though from a purely formal view- 
point it is indeed possible that a dimensionless coefficient 
like C x can be any function of Ri gra d ■ 

To expose the physical reason for the failure of the 
down-gradient approach, recall that in a stratified flow, 
in the presence of gravity, the turbulent kinetic energy 
is not an integral of motion. Only the total mechanical 
energy, the sum of the kinetic and the potential energy, is 
conserved in the inviscid limit. As it was shown already 
by Richardson, the difficult point is that an important 
contribution to the potential energy comes not just from 
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the mean density profile, but from the density fluctu- 
ations. Clearly, any reasonable model of the turbulent 
boundary layer must obey the conservation laws. 

The physical requirement of conserving the total me- 
chanical energy calls for an explicit consideration not 
only of the mean profiles, but also of all the relevant 
second-order, one-point, simultaneous correlation func- 
tions of all the fluctuating fields together with some clo- 
sure procedure for the appearing third order moments. 
First of all, in order to account for the important ef- 
fect of stratification on the anisotropy, we must write 
explicit equations for the entire Reynolds stress tensor, 
Tij = (uiUj) . Next, in the case of the temperature strat- 
ified turbulent boundary layers we follow tradition [see, 
e.g. Zeeman (1981), Hunt et al. (1988), Schumann and 
Gerz (1995), Hanazaki and Hunt (2004), Keller and van 
Atta (2000), Stretch et al. (2001), Elperin ct al. (2002), 
Cheng et al. (2002) Luytcn ct al. (2002), and Rchmann 
and Hwang (2005)] and account for the turbulent poten- 
tial energy which is proportional to the variance of the 
potential temperature deviation, (9 2 ). And last but not 
least, we have to consider explicitly equations for the ver- 
tical fluxes, t xz and F z , which include the down gradient 
terms proportional to the velocity and temperature gra- 
dients, and counter-gradient terms, proportional to F x 
(in the equation for t xz ) and to (# 2 ) (in the equation for 
F z ). 

Unfortunately, the resulting second order closure seems 
to be inconsistent with the variety of boundary-layer 
data, and many authors took the liberty to introduce 
additional fitting parameters and sometimes fitting func- 
tions to achieve a better agreement with the data (see re- 
views of Umlauf and Burchard (2005), Weng and Taylor 
(2003), Zeeman, (1981), Melor and Yamada (1974), and 
references therein). Moreover, in the second order clo- 
sures the problem of critical Richardson number seems 
to persists (Cheng et al., 2002; Canuto, 2002). 

Notice that in spite of obvious inconsistency of the 
first-order schemes, most of the practically used tur- 
bulent models are based on the concept of the down- 
gradient transport. One of the reasons is that in the 
second-order schemes instead of two down-gradient equa- 
tions (0.3) one needs to take into account eight nonlinear 
coupled additional equations i.e. four equations for the 
Reynolds stresses, three equations for the heat fluxes and 
equation for the temperature variance. As the result, the 
second-order schemes have seemed to be rather cumber- 
some for comprehensive analytical treatment and have 
allowed to find only some relationships between correla- 
tion functions (see, e.g., Cheng et al., 2002). Unfortu- 
nately, the numerical solutions to the complete set of the 
second-order schemes equations which involve too many 
fitting parameters are much less informative in clarifica- 
tion of physical picture of the phenomenon than desired 
analytical ones. 

In this paper we suggest a relatively simple second- 
order closure model of turbulent boundary layer with 
stable temperature stratification that, from one hand, 



accounts for main relevant physics in the stratified TBL 
and, from the other hand, is simple enough to allow com- 
plete analytical treatment including the problem of crit- 
ical Rigrad- To reach this goal we approximate the third 
order correlations via the first- and second-order ones, ac- 
counting only for the most physically important terms. 
We will try to expose the approximations in a clear and 
logical way, providing the physical justification as we go 
along. Resulting second-order model consist of nine cou- 
pled equations for the mean velocity and temperature 
gradients, four components of the Reynolds stresses, two 
components of the temperature fluxes and the tempera- 
ture variance. Thanks to the achieved simplicity of the 
model we found an approximate analytical solution of 
these equations, expressing all nine correlations as func- 
tions of only one governing parameter, £{z)/L 1 where £(z) 
is the outer scale of turbulence (depending on the eleva- 
tion z and also known as the "dissipation scale") and L 
- is the Obukhov length. 

We would like also to stress, that in our approach 
£(z)/L is an external parameter of the problem. For 
small elevations z -C L, it is well accepted that £{z) is 
proportional to z, while the £(z) dependence is still under 
debate for z comparable or exceeding L. For z > L the 
assignment and discussion of the actual dependence of 
the outer scale of turbulence, £(z), which is manifested 
in the nature is out of the scope of this paper, and is 
remained for future work. At time being, we can analyze 
consequences of our approach for the following versions 
of £{z) dependence at z » L: 

• function £(z) is saturated at some level of the order of 
L. For concreteness we can take 

l/t(z) = V(diz)- 2 + (d 2 L)- 2 , di~d2~l. (0.4) 

• £(z) is again proportional to z for elevations much 
larger than L: £{z) — d 3 z but with the proportionality 
constant d 3 < d\ . If so, we can also study the case £(z) ^> 
L even though such a condition may not be realizable in 
Nature. In that case our analysis of the limit £(z) 3> L 
has only a methodological character: it allows to derive 
an approximate analytic solution for all the objects of 
interest as functions of £{z)/L that is also valid for the 
outer scale of turbulence not exceeding a value of the 
order of L. 

It should be noticed that traditional turbulent closures 
(including ours) cannot be applied for strongly strati- 
fied flows with Rigrad ^ 1 (may be even at Ri gra d ~ 1). 
The main reason is that these closures are roughly jus- 
tified for developed vortical turbulence, in which the 
eddy-turnover time is of the order of its life time; in 
other words, there are no well defined "quasi-particlcs" 
or waves. This is not the case for stable stratification 
with Ri g rad ^ 1, in which the Brunt- Vaisala, frequency 

N= y/0de(z)/dz, (0.5) 

is larger then the eddy-turnover frequency 7. It means 
that for Rigrad }Z 1 there are weakly decaying Kelvin- 
Helmoholtz internal gravity waves (with characteristic 
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frequency iV and decay time above 1 /j) , propagating on 
large distances, essentially effecting on TBL, as pointed 
out by Zilitinkevich, (2002). We concentrate in our pa- 
per on self-consistent description of the lower part of the 
atmospheric TBL, in which turbulence has vortical char- 
acter and consequently, large values of Ri gra d do not ap- 
pear. We relate large values of Ri gra d in the upper part 
of TBL with contributions of the internal gravity waves 
to the energy and the energy flux in TBL, to the momen- 
tum flux, and to the production of (vortical) turbulent 
energy. Due to their instability in a shear flow, the waves 
can break and create turbulent kinetic energy. All these 
effects are beyond the framework of our paper. Their 
description in the upper "potential-wave" TBL and in- 
termediate region with the combined "vortical-potential" 
turbulent velocity field is in our nearest agenda. 

To make the paper more transparent for wide audience, 
not necessarily experts in atmospheric TBL, we attempt 
to present the material in a self-contained manner, and 
organized it as follows. 

In Sect. IA we use the Oberbeck-Boussinesq approxi- 
mation and apply the standard Reynolds decomposition 
(into mean values and turbulent zero-mean fluctuations 
of the velocity and temperature fields) to derive equa- 
tions for the mean values and balance equations for all 
relevant second-order correlation functions. In Sect. IB 
we demonstrate that the resulting balance equations ex- 
actly preserve (in the non-dissipative limit) the total me- 
chanical energy of the system, that consists of the kinetic 
energy of the mean flow, kinetic energy and potential en- 
ergy of the turbulent subsystem. 

In Sect. II we describe the proposed closure proce- 
dure that results in a model of stably stratified TBL, 
that accounts explicitly for all relevant second-order cor- 
relations. The third order correlations which appear in 
the theory are modeled in terms of second-order correla- 
tions in Sects. IIA and B. Further simplifications are pre- 
sented in Sects. IIC and D for stationary turbulent flows 
in a plane geometry outside the viscous and buffer layers. 
In Sect. HE we suggest a generalization of the standard 
"wall-normalization" to obtain the model equations in a 
dimensionless form with only one governing parameter, 
£{z)/L. 

Section III A contains approximate analytical solution 
of the model. It is shown that the analytical solution 
deviates from the numerical counterpart in less than a 
few percent in the entire interval < {£/ L) < oo. 

The last Sect. Ill is devoted to a detailed description 
of our results: profiles of the mean velocity and potential 
temperature (Sect. IIIA), profiles of the turbulent kinetic 
and "temperature" energies, profiles of the anisotropy of 
partial kinetic energies (Sect. IIIB), profiles of the tur- 
bulent transport parameters v T and % T , profiles of the 
gradient- and flux-Richardson numbers Ri gra d and Rifl ux , 
and the dependence of the turbulent Prandtl number Pr T 
vs. i/L and Ri gra d, Sect. IIIC. In conclusion Sect. HID, 
we consider the validity of the down-gradient transport 
concept (0.3) and explain why it is violated in the up- 



per part of TBL. The problem of critical Ri g r a d is 

also 

discussed. 



I. SIMPLIFIED DYNAMICS IN A STABLY 
TEMPERATURE-STRATIFIED TBL AND THEIR 
CONSERVATION LAWS 

The aim of this section is to consider the simplified 
dynamics of a stably temperature-stratified turbulent 
boundary layer, aiming finally at an explicit description 
of the height dependence of important quantities like the 
mean velocity, mean temperature, turbulent kinetic and 
potential energies, etc. In general one expects very differ- 
ent profiles from those known in standard (unstratified) 
wall-bounded turbulence. We want to focus on these dif- 
ferences and propose that they occur already relatively 
close to the ground allowing us to neglect (to the leading 
order) the dependence of the density on height and the 
Coriolis force. We thus begin by simplifying the hydro- 
dynamic equations which are used in this section. 



A. Simplified hydrodynamic equations and 
Reynolds decomposition 

First we briefly overview the derivation of the gov- 
erning equations in the Boussinesq approximation. The 
system of hydrodynamic equations describing a fluid in 
which the temperature is not uniform consists of the 
Navier-Stokes equations for the fluid velocity, U(r,t), 
a continuity equation for the space and time depen- 
dent (total) density of the fluid, p(r,t), and of the heat 
balance equation for the (total) entropy per unit mass, 
S(r,t), Landau and Lifshitz, 1987. 

These equations are considered with boundary condi- 
tions that maintain the solution far from the equilibrium 
state, where U = S = 0. These boundary conditions are 
U = at zero elevation, U = const at a high elevation of 
a few kilometers. This reflects the existence of a wind at 
high elevation, but we do not attempt to model the phys- 
ical origin of this wind in any detail. The only important 
condition with regards to this wind is that it maintains a 
momentum flux towards the ground that is prescribed as 
a function of the elevation. Similarly, we assume that a 
stable temperature stratification is maintained such that 
the heat flux towards the ground is prescribed as well. 

We neglect the viscous entropy production term as- 
suming that the temperature gradients are large enough 
such that the thermal entropy production term domi- 
nates. For simplicity of the presentation we restrict our- 
selves by relatively small elevations and disregard the 
Coriolis force (for more details, see Wyngaard, 1992). 
On the other hand we assume that the temperature and 
density gradients in the entire turbulent boundary layer 
are sufficiently small to allow employment of local ther- 
modynamic equilibrium. In other words, we assume the 
validity of the equation of state. 
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As a "basic reference state" (BRS) denoted hereafter 
by a subscript 'V we use the isentropic model of the 
atmosphere, where the entropy is considered space ho- 
mogeneous. Now assuming smallness of deviations of the 
density and pressure from their BRS values and exploit- 
ing the equation of state, one obtains a simplified equa- 
tion, which is already very close to the standard Navier- 
Stokes equation in the Boussinesq approximation. Intro- 
ducing (generalized) potential temperature, one results in 
the well-known system of hydrodynamic equations in the 
Boussinesq approximation. Close to the ground, where 
one can neglect the dependence of the density on height, 
the system reads: 

f = -^-/30 d + ,A W ,^ = x Ae d .(l.l) 
Ut p\> ut 

Here V/Vt = d/dt + U- V is the convection time deriva- 
tive, p - deviation of pressure from BRS, pb is the density 
in BRS, (3 = gfi is the buoyancy parameter (f3 = — z(3, 
[3 = g/3, g is the gravity acceleration and (3 is the ther- 
mal expansion coefficient, which is equal to 1/T, recip- 
rocal molecular temperature, for an ideal gas), 6d is the 
deviation of the potential temperature from BRS value, 
v - kinematic viscosity and x is the kinematic thermal 
conductivity. 

To develop equations for the mean quantities and cor- 
relation functions one applies the Reynolds decompo- 
sition: u = u + u, (u) = u, (u) = o , e d = 
Q + e, (e d ) = e, (e) = o, P = ( P ) +p, (p) = o. 

Here the average (• • • ) stands for an averaging over a 
horizontal plane at a constant elevation. This leaves the 
average quantities with a z,t dependence only. Substi- 
tuting in Eqs. (1.1) one gets equations of motion for the 
mean velocity and mean temperature profiles 



Dt 



Pb 



Dt 



Here D/ Dt = d/dt+U -V is the mean convection deriva- 
tive. The total (molecular and turbulent) momentum 
and thermal fluxes are 



-vVj Ui 



1 1] , 



f = -xve 



(1.3) 



where = (iiiUj) is the Reynolds stress tensor describ- 
ing the turbulent momentum flux, and F = (u8) is the 
turbulent thermal flux. In order to derive equations for 
these correlation functions, one considers the equations 
of motion for the fluctuating velocity and temperature: 

Du/Dt = — u • VU -u-Vu+(u- Vu) (1.4a) 

-(Vp/p b ) + vAu- (3 9, 
D8/Dt = -u- V9-U- V9 + X A0 + (u- V0).(1.4b) 

The whole set of the second order correlation functions 
includes the Reynolds stress, r^ , the turbulent thermal 
flux, F, and the "temperature energy" Eg = (# 2 ) /2, 
which is denoted and named by analogy with the tur- 
bulent kinetic energy density (per unit mass and unit 



volume), E K = (\u\ 2 )/2 = Tr{ Tij }/2. Using (1.4) one 
gets the following "balance equations" : 

jj ^ &ij ~t~ ^ijk + Rij , (1.5a) 



DF\ 
Dt 
DE e 
Dt 



d 

"i" "7^ y\i + Si , 

OXi 



+ e + V T 



-F ■ V6 . 



(1.5b) 
(1.5c) 



Here we denoted the dissipations of the Reynolds-stress, 
heat-flux and the temperature energy by 



e 



\dx k dx k I 
X<|V0| 2 ), 



\ dx k dx k I ' 



(1.6) 



The last term on the LHS of each of Eqs. (1.5) describes 
spatial flux of the corresponding quantity. In models of 
wall bounded unstratified turbulence it is known that 
these terms are very small almost everywhere. We do not 
have sufficient experience with the stratified counterpart 
to be able to assert that the same is true here. Neverthe- 
less, for simplicity we are going to neglect these terms. 
It is possible to show that the accounting for these terms 
does not influence much the results. Note that keep- 
ing these terms turns the model into a set of differential 
equations which are very cumbersome to analyze. This 
is a serious uncontrolled step in our development, so we 
cross our fingers and proceed with caution. Since these 
terms are neglected we do not provide here the explicit 
expressions for T^, Tjj, and T. 

The first term on the RHS of the balance Eq. (1.5a) for 
the Reynolds stresses is the "Energy Production tensor" 
Vij, describing the production of the turbulent kinetic 
energy from the kinetic energy of the mean flow, propor- 
tional to the gradient of the mean velocity: 



Vij = -Tj fe dUj/dxk - T jk dlli/dxk ■ 



(1.7a) 



The second term on the RHS of Eq. (1.5a), Cij, will be 
referred hereafter to as the "Energy Conversion tensor" . 
It describes the conversion of the turbulent kinetic energy 
into potential energy. This term is proportional to the 
buoyancy parameter (3 and the turbulent thermal flux F: 



dj^-ptFiSjz + FjSiz) . 



(1.7b) 



The next term in the RHS of Eq. (1.5a) is known as the 
"Pressure-rate-of-strain tensor" : 

TZij = (ps i:j /p h ) , Sij = dui/dxj + duj/dxi . (1.7c) 

In incompressible turbulence its trace vanishes, therefore 
TZij does not contribute to the balance of the kinetic en- 
ergy. As we will show in Sec. II A, this tensor can be pre- 
sented as the sum of three contributions (Zeman, 1981), 



dRI i p ip _j_ pi 



(1.7d) 
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in which Rf? is responsible for the nonlinear process of 
isotropization of turbulence and is traditionally called the 
"Return-to-Isotropy" , R-J is similar to the energy pro- 
duction tensor (1.7a) and is called "Isotropization of Pro- 
duction". A new term, appearing in the stratified flow, 
R-j, is similar to the energy conversion tensor (1.7b) and 
will be refereed to as the "Isotropization of Conversion" . 

Consider the balance of the turbulent thermal flux F, 
Eq. (1.5b). The first term in the RHS, A., describes the 
source of F and, by analogy with the energy-production 
tensor, Vij, is called "Thermal-flux production vector". 

Like Vij, Eq. (1.7b), it has the contribution, A i , pro- 
portional to the mean velocity gradient: 



A = A l +4 +A 



E6 
i > 



A, = -F-VUi 



The balance Eq. for E K follows from Eq. (1.2): 

DE K /D t + v (VjUrf + V 3 {Ui nj) = 

[source E/c] + TijVj Ui , (1.8a) 

with the help of identity: UiVj Tij = Vj (UiTij) — TijVj Ui 
and definition (1.3). The terms on the LHS of this Eq., 
proportional to v and respectively, describe the dissi- 
pation and the spatial flux of Ejc- The term [source E/c] 
on the RHS of Eq. (1.8a) describes the external source 
of energy, originating from the boundary conditions de- 
scribed above, and Tij'Vj Ui describes the kinetic energy 
out-flux from the mean flow to turbulent subsystem. 

The balance Eq. for the turbulent kinetic energy fol- 
lows directly from Eq. (1.5a): 



-TijdS/dxj, A, = 2j3E e 5 lz , (1.7e) DE K /Dt + [s u + Vj T Uj ]/2 = -t^ Ui + (3F Z . (1.8b) 



and two additional contributions, related to the temper- 
ature gradient and to the "temperature energy", Eg, and 
the buoyancy parameter. One sees, that in contrary to 
the oversimplified assumption (0.3b), the thermal flux in 
such a turbulent media cannot be considered as propor- 
tional to the temperature gradient. It has also a con- 
tribution proportional to the velocity gradient and even 
to the square of the temperature fluctuations. Moreover, 
the RHS of the flux-balance Eq. (1.5b) has an additional 
term, the "Pressure-temperature-gradient vector" which, 
similarly to the pressure-rate-of-strain tensor (1.7d), can 
be divided into three parts (Zeman, 1981): 



B = (pVe/p b ) = B RD + B sv + B E 



(1.7f) 



As we will show in Sec. II A the first contribution, Bf° cx 
(u u Vj 9) is responsible for the nonlinear flux of F in 
the space of scales toward smaller scales, similarly to the 
correlation (uuu), which is responsible for the flux of 
kinetic energy (u 2 }/2 toward smaller scales. The corre- 
lation Bf ° cx (uu\7i 9) may be understood as the non- 
linear contribution to the dissipation of the thermal flux. 
Correspondingly we will call it "Renormalization of the 
Thermal-Flux Dissipation" and will supply it with a su- 
perscript " RD " . The next two terms in the decomposi- 
tion (1.7f) are B i cx S v and B i cx Eg. They describe 
the renormalization of the thermal-flux production terms 
A i cx and A i cx Eg, accordingly. 



B. Conservation of total mechanical energy in the 
exact balance equations 

The total mechanical energy of temperature stratified 
turbulent flows consists of three parts with densities (per 
unit mass): E = Ejc+E k + E p , where E K = \U\ 2 /2 is the 
density of kinetic energy of the mean flow, E K = tu/2 is 
the density of turbulent kinetic energy and E P — f3Eg / Sq 
is the density of potential energy, associated with turbu- 
lent density fluctuation p = (3 9 pi, , caused by the (poten- 
tial) temperature fluctuations 9, and S® = dO/dz. 



On the LHS of Eq. (1.8b) one sees the dissipation and 
spatial flux terms. The first term on the RHS originates 
from the energy production, | Va, defined by Eq. (1.7a). 
This term has an opposite sign to the last term on the 
RHS of Eq. (1.8a) and describes the production of the 
turbulent kinetic energy from the kinetic energy of the 
mean flow. The last term on the RHS of Eq. (1.8b) origi- 
nates from the energy conversion tensor \ Ca, Eq. (1.7b), 
and describes the conversion of the turbulent kinetic en- 
ergy into potential one. 

According to the last of Eqs. (1.5), one gets the bal- 
ance equation for the potential energy E P ; multiplying 
Eq. (1.5c) for Eg by (3/S @ : 



DE P /Dt + (3 e + VjTj /S@ = -f3F : 



(1.8c) 



The RHS of this Eq. [coinciding up to a sign with the last 
term on the RHS of Eq. (1.8b)] is the source of potential 
energy (from the kinetic one). 

In the sum of the three balance equations, the conver- 
sion terms (of the kinetic energy from the mean to tur- 
bulent flows and of the turbulent kinetic energy to the 
potential one) cancel and one gets the total mechanical 
energy balance: 

DE /Dt+ [diss E] + V [flux£] = [source E K ] . (1.9) 

This equation exactly respects the conservation of total 
mechanical energy in the dissipation- less limit, irrespec- 
tive of the closure approximations. This is because the 
energy production and conversion terms are exact and do 
not require any closures, while the pressure-rate-of-strain 
tensor, that requires some closure, does not contribute to 
the total energy balance. 



II. THE CLOSURE PROCEDURE AND THE 
RESULTING MODEL 

In this section we describe the proposed closure proce- 
dure that results in a model of stably stratified TBL. In 
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developing this model we strongly rely on the analogous 
well developed modeling of standard (unstratified) TBL. 
The final justification of this approach can be done only 
in comparison to data from experiments and simulations. 
We will do below what we can to use the existing data, 
but we propose at this point that much more experimen- 
tal and simulational work is necessary to solidify all the 
steps taken in this section. 



A. Pressure-Rate-of-Strain tensor IZij and 
Pressure-Temperature-Gradient vector B 

The correlation functions IZij and 0, de- 
fined by Eqs. (l-7c) and (1.7f), include fluc- 
tuating part of the pressure p. The Pois- 
son's equation for p follows from Eq. (1.4): 



Ap = p h 



-VjVj (uiUj - (uiUj) + UiUj + UjUi)+(3V z 9 



The solution of this equation includes a harmonic part, 
Ap = 0, which is responsible for sound propagation 
and does not contribute to turbulent dynamics at 
small Mach numbers. Thus this contribution can be 
neglected, the inhomogeneous solution includes three 
parts p = p b [Puu +PUu +Pe], where 

p uu = A _1 VjVj ((uiUj) - UiUj) , (2.1) 
p Uu = A^ViVj {U iUj + U jUi ) , Pe = PA^V Z 9 , 

and the inverse Laplace operator A -1 is defined as usual 
in terms of an integral over the Green's function. 

Correspondingly the correlations IZij and B consist of 
three terms, Eqs. (1.7d) and ( 1 . 7f ) , in which 

Rfj = (PuuSij) , Rf, = (pun Sij) , Rfj = (pesij) ,(2.2) 



Bf 



{PuuVe) , b = {p Uu ve) , B i = ( Pe ve) . 



All of those terms originating from p uu are the most prob- 
lematic because they introduce coupling to triple corre- 
lation functions: Rfj cx {uiUjUk) and B RD cx (u 2 W9). 
Thus they require closure procedures whose justification 
can be only tested a-posteriori against the data. 

Having in mind to simplify the model in most possible 
manner, we adopt for the diagonal part of the Return- 
to-Isotropy tensor, the simplest Rota form (Rotta, 1951) 



(2.3a) 



in which 7 RI is the relaxation frequency of diagonal com- 
ponents of the Reynolds-stress tensor toward its isotropic 
form, 2E K /3. The parametrization of 7 M will be dis- 



with, generally speaking, 7 M ^ 7 RI . Moreover, on the in- 
tuitive level, we can expect that off-diagonal terms should 
decay faster then the diagonal ones, i.e. 7 RI > 7 RI . In- 
deed, our analysis of DNS results shows that 7 w /7w — 
1.46 (L'vov et al., 2006b). 

The term B RO also describes return-to-isotropy due to 
nonlinear turbulence self interactions (Zeman, 1981), and 
may be modeled as: 



Bf 



(2.3c) 



This equation dictates the vectorial structure of Bf n oc 
Fi, which will be confirmed below. The rest can be un- 
derstood as the definition of the 7rd as the relaxation 
frequency of the thermal flux. Its parametrization is the 
subject of further discussion in Sec. II D. 

The traceless "Isotropization-of-Production" tensor, 
RIJ, has a very similar structure to the production ten- 
sor, Vij, Eq. (1.7a), and thus is traditionally modeled in 
terms of Vij (Pope, 2001): 

R% ~ -C IP (Pij - % P/3) , P = Ti {Pij } • (2.3d) 

The accepted value of the numerical constant C IP = | 
(Pope, 2001). 

The traceless "Isotropization-of- Conversion" tensor, 
RJf does not exist in unstratified TBL. Its structure is 
very similar to the conversion tensor, Cij, Eq. (1.7b). 
Therefore it is reasonable to model it in the same way 
in terms of (Zeman, 1981): 



R 



„ ~ -C m (Cij - % C/3) , C = Tr {dj} , (2.3e) 



with some new constant Gc ■ 

The renormalization of production terms B i and B i 
are very similar to the corresponding thermal flux pro- 
duction terms, A i and A i , defined by Eqs. (1.7e). 
Therefore, in the spirit of Eqs. (2.3d) and (2.3c), they 
are modeled as follows: 

B. U = (C su - 1)4" - (1 - C SU )(F -V)Ui, (2.3f) 

EO EQ 

Bi = -(C Be +l)Ai =-2p(C B0 +l)E e 5 iz .(2.3g) 
Using this and (2.1) one finds the sign of C E0 : 

-p(C Be +l)E e = {peV z 6) = /3({V Z 6)A-\V Z 6)), 

c Ee = (i + {{v z e)A-\v z e))/{e 2 )) < o . (2.4) 



To estimate C Eg we assume that on the gradient scales 
the temperature fluctuations are roughly isotropic, and 
therefore we can estimate A = V 2 X + + V 2 Z w 3V^. 
Introducing this estimate and integrating by parts leads 



cussed later. The tensor R-V is traceless, therefore the to C E 
frequency 7 RI must be the same for all the diagonal com- 
ponents of On the other hand there are no reasons 
to assume that off-diagonal terms have the same relax- 
ation frequency. Therefore, following L'vov et al. (2006a) 
we assume that 



-2/3. 



Reynolds-stress-, thermal-flux-, and 
thermal-dissipation 



R™ ~ -7 RI r i: 



Far away from the wall and for large Reynolds num- 
(2.3b) bers the dissipation tensors are dominated by the viscous 
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scale motions, at which turbulence can be considered as 
isotropic. Therefore, the vector e should vanish, while 
the tensor e^, Eq. (1.6), should be diagonal: 

ej = 0, Sij =2'y uu E K 6 ij /3, (2.5a) 

where the numerical prefactor | is chosen such that j uu 
becomes the relaxation frequency of the turbulent kinetic 
energy. Under stationary conditions the rate of turbu- 
lent kinetic energy dissipation is equal to the energy flux 
through scales, that can be estimated as {uuu} /£, where 
I is the outer scale of turbulence. Therefore, the natu- 
ral estimate of 7 Mtt involves the triple- velocity correlator, 
luu ~ ((uuu) /£ (uu)), exactly in the same manner, as the 
Return-to-Isotropy frequencies, 7 RI and7 RI inEqs. (2.3a) 
and (2.3b). Similarly, 

e = l00 E e , lee ~ (98u) /£ (69) . (2.5b) 



C. Stationary balance equations in plain geometry 

In the plane geometry, the equations simplify fur- 
ther. The mean velocity is oriented in the (stream- 
wise) x direction and all mean values depend on the 
vertical (wall-normal) coordinate z only: U = U(z) x, 

9 = Q(z), Ti j = Tij(z), F = F(z), Eg = Eg(z). 

Therefore (U ■ V) (. . . ) =0, and in the stationary case, 
when d jdt = 0, the mean convective derivative van- 
ishes: D /Dt = 0. Moreover due to the y — > — y sym- 
metry of the problem the following correlations vanish: 
T xy = T yz — Fy = 0. The only non-zero components of 
the mean velocity and temperature gradients are: 

S u = dU/dz , S B = dQ/dz . (2.6) 



1. Equations for the mean velocity and temperature profiles 

Having in mind Eqs. of Sec. II C and integrating Eqs. 
(1.2) for U x and O over z, one gets equations for the 
total (turbulent and molecular) mechanical-momentum 
flux, t(z), and thermal flux, F, toward the wall 

t xz (z) = -vS v + t xz T xz (0) = -r* , (2.7a) 
F z (z) = - x S e +F z ^ F z (0) = -F, . (2.7b) 

The total flux of the x-component of the mechanical mo- 
ment in z-direction is Pbf xz (z) = J dz(d (p) jdx) + const. 
Generally speaking, t xz (z) depends on z. For example, 
for the pressure driven planar channel flow (of the half- 
wight S) p b r xz (z) = (d (p) /dx)(S -z)<0. 

Relatively close to the ground, where z <C 6, the z de- 
pendence of t xz (z) can be neglected. In the absence of 
the mean horizontal pressure drop and spatial distributed 
heat sources r and F are z-indepcndent, and thus equal 
to their values at zero elevation, as indicated in Eqs. (2.7) 



after "=>"-sign. Notice, that in our case of stable strati- 
fication both vertical fluxes, the a;-component of the me- 
chanical momentum, t XZi and the thermal flux, F z , are 
directed toward the ground, i.e. negative. For the sake 
of convenience, we introduce in Eqs. (2.7) notations for 
their (positive) zero level absolute value: r* and . 

Recall that in the plain geometry U z — 0. Nevertheless 
one can write an equation for U z : 

d(T zz + (p)/p b )/dz = (3Q, (2.8) 

which describes a turbulent correction (cx r zz ) to the hy- 
drostatic equilibrium. Actually, this equation determines 
the profile of (p) , that does not appear in the system of 
balance equations (2.7). 

2. Equations for the pair (cross)-correlation functions 

Consider first the balance Eqs. (1.5a) for the diago- 
nal components of the Reynolds-stress tensor in algebraic 
model (which arises when we neglect the spatial fluxes) : 

TE K + 3 7RI r ra /2 = - (3 - 2 C IP ) t xz S v - C IC f3 F z , 
TE K + 3 7RI r w /2 = -C ip t xz S v - C 1C 13 F z , (2.9) 
TE K + 3j k1 t zz /2 = —C' ip t xz S u + 

where F = "f uu — 7 RI . The LHS of these equations 
includes the dissipation and Return-to-isotropy terms. 
On the RHS we have the kinetic energy production and 
isotropization of production terms (both proportional to 
S u ) together with the conversion and isotropization of 
conversion terms, that are proportional to the vertical 
thermal flux F z . The horizontal component of the ther- 
mal flux, F x , does not appear in these equations. 

System (2.9) allows to find anisotropy of the turbulent- 
velocity fluctuations and to get the balance Eqs. for the 
turbulent kinetic energy with the energy production and 
conversion terms on the RHS: 

3t xx = 2{[2(l-a P )W7Ri + l]^K (2.10a) 

-(3-2C 1P + C lc )(3F z / 7m }, 
3T yy = 2{[(C I p-l)r u „/7 RI + l]£ K (2.10b) 

-(c IP + c IC )/3iy 7RI }, 

3t zz = 2{[(C IP -l)r iiU /7 M + l]£ K (2.10c) 
-(C IP -2C ic -3)/3F z / 7ri }, 
T UU E K = -T xz S L ,+f3F z , (2.10d) 

Equation (2.10d) includes the only non- vanishing tan- 
gential (off-diagonal) Reynolds stress t xz and has to be 
accompanied with an equation for this object: 

7mr„ = (C IP - l)r„ S v + (1 + C lc )f3F x . (2.10e) 

This equation manifests that the tangential Reynolds 
stress t xz , that determines the energy production [ac- 
cording to Eq. (2.10d)], influences, in its turn, on the 
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value of the streamwise thermal flux F x , which therefore 
effects on the turbulent kinetic energy production. 

As we mentioned, in the plain geometry F y = 0. Equa- 
tions (1.5b) for the F x and F z in this case take the form: 

i™F x = -(t xz S + C sv F z S u ) , (2.11a) 

lnnF z = -(T zz S e + 2C Ee /3Eg) , (2.11b) 

in which the RHS describes the thermal-flux production, 
corrected by the isotropization of production terms. 

The last Eq. (1.5c) for Eg, represents the balance be- 
tween the dissipation (LHS) and production (RHS): 

lee E e = -F z S e . (2.11c) 



D. Simple closure of time-scales and the balance 
equations in the turbulent region 

At this point we follow a tradition in modeling of all the 
nonlinear inverse time-scales by dimensional estimates 
(Kolmogorov, 1941): 

luu = Cuuy^/e, 7 RI = C M 7„ U , (2.12) 
7w = Cw7w j lee = Cggj uu , 7rd = C u eluu ■ 

Remember that I is the "outer scale of turbulence" . This 
scale equals to z for z < L, where L is the Obukhov 
length (definition is found below). 

Detailed analysis of experimental, DNS and LES data 
(see L'vov et al., 2006, and references therein) shows that 
for unstratified flows, g = 0, the anisotropic boundary 
layers exhibits values of the Reynolds stress tensor that 
can be well approximated by the values t xx — E K , t vv = 
t zz = E K /2. In our approach this dictates the choice 
C RI = 4(1 — C IP ). Also we can expect that r yy is almost 
not affected by buoyancy. This gives simply C IC = — C IP . 
If so, Eqs. (2.10) with the parametrization (2.12) can be 
identically rewritten as follows: 

TxX ~ K 2 7M n' TyV ~ 2 ' T "~ 2 + 2 7 ™' 
JuuE K = f3F z - t xz S v , 7„ u = c uu \[~E^I I , (2.13a) 
4 C M -f uu T xz = (3 F x - t zz S v . 

For completeness we also repeated here the parametriza- 
tion (2.12) of 7„„. Finally we present the version of the 
balance Eqs. for the thermal flux (2.11a), (2.11b), and 
for the "temperature energy", (2.11c), after all the sim- 
plified assumptions: 

CeeluuEg = —F z S e , 

C u9 luuF x = - (r xz S e + C SU F Z S a ) , (2.13b) 
C u eluuF z — — {r zz S e + 2C Ee 0Eg) . 



E. Generalized wall normalization 

The analysis of the balance Eqs. (2.13) is drastically 
simplified if they are presented in a dimcnsionless form. 
Traditionally, the conventional "wall units" are intro- 
duced via the wall friction velocity = y 7 ?*, and the 
viscous length-scale A„ = v/u*. A wall unit for the tem- 
perature 9* = F*/u* is defined via the thermal flux at 
the wall and friction velocity. Subsequently, r + = r/A*, 
t+ = tK/u„ U+ = U/u*,p+ = p/p b ul 9+ = Q/0„ 
6 + = 9/9*, etc. Then the governing Eqs. (1.1) take the 
form: 

V+U + /Vt+ + V+p+ = zO+/L+ +A+U+, 

V+e+/Vt + = A+6+/Pr. (2.14) 

These Eqs. include two dimensionless parameters: the 
conventional Prandtl number Pr= vjn, and L + the 
Obukhov length L measured in wall units: L = u^/f3F*, 
L + = L/A*. Wc used here the modern definition of the 
Obukhov length, which differs from the old one by the 
absence of the von-Karman constant k in its denominator 
(Monin and Obukhov, 1954). 

Outside of the viscous sub-layer, where the kinematic 
viscosity and kinematic thermal conductivity can be ig- 
nored, L + is the only dimensionless parameter in the 
problem, which separates the region of weak stratifica- 
tion, z + < L + , and the region of strong stratification, 
where z + > L + . 

Given the generalized wall normalization we introduce 
objects with a superscript " + " in the usual manner: 

S+ = US V , S+ = K S e /9* , 7+ = ^7 , t± = nj/ul , 
F+ = F/u*6* , E+ = Eg /6l . (2.15) 

In the turbulent region, governed by L + only, Eqs. (2.7) 
simplify to r+, = -1, F+ = -1. 



F. Rescaling symmetry and ^-representation 

Outside of the viscous region, where Eqs. (2.13) were 
derived, the problem has only one characteristic length, 
i.e. the Obukhov scale L. Correspondingly, one expects 
that the only dimensionless parameter that governs the 
turbulent statistics in this region should be the ratio of 
the outer scale of turbulence, i, to the Obukhov length- 
scale L, which we denote as ft = l/L = C + / L + . Indeed, 
introducing "|-objects" : 

£i = l/L , S\ = S+1+ , Si = S+£+ , (2.16) 
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and using Eqs. (2.15) one rewrites the balance Eqs. (2.13) 
as follows: 



T XX ~ + ^ /2C UU 






(2.17a) 


2r+=E+-£*/c uu 






(2.17b) 


r R + 1 - ftF + - 


' xz^u 5 




(2.17c) 


4 C RI c uu \J EkT^z = 


- 


T ZZ D u 7 


(2.17d) 


Cee c U u \J E^Eg = - 






(2.17c) 


C u6 c uu sJe+F+ = - 




-CsvF+S*, 


(2.17f) 


Cuec uu \[~E+F+ = - 




— 2 C Ee ft E^ , 


(2-17g) 



These equations are the main result of current Sec. II. It 
may be considered as "Minimal Model" for stably strat- 
ified TBL, that respects the conservation of energy de- 
scribes anisotropy of turbulence and all relevant fluxes 
explicitly and, nevertheless is still simple enough to allow 
comprehensive analytical analysis, that results in an ap- 
proximate analytical solution (with reasonable accuracy) 
for the mean velocity and temperature gradients S v and 

5 , and all second-order (cross)-correlation functions. 
As expected, the only parameter appearing in the Min- 
imal Model (2.17) is ft. The outer scale of turbulence, 
£, does not appear by itself, only via the definition of ft 
(2.16). Therefore our goal now is to solve Eqs. (2.17) 
in order to find five functions of only one argument ft: 

51, Si, E+, E+ and F+ . After that we can specify 
the dependence C + (z + ) and then reconstruct the z + - 
dependence of these five objects. 



III. RESULTS AND DISCUSSION 

A. Analytical solution of the Minimal-Model 
balance equations (2.17) 

This subsection is devoted to an analytical and numer- 
ical analysis of the Minimal- Model (2.17). An example of 
numerical solution of Eqs. (2.17) (with some reasonable 
choice of the phenomenological parameters) is shown in 
Fig. 1. Nevertheless, it would be much more instructive 
to have approximate analytical solutions for all corre- 
lations that will describe their ^-dependence with rea- 
sonable accuracy. The detailed cumbersome procedure 
of finding these solutions is skipped here, but a brief 
overview is as follows. 

The Eqs. (2.17) can be reformulated as a polynomial 
equation of ninth order for the only unknown \J E£ ■ An 
analysis of its structure helps to formulate an effective 
interpolation formula (3.1), discussed below. Hence, we 
found the solutions of Eqs. (2.17) at neutral stratifica- 
tion, ft = 0, corrected up to the linear order in ft. Its 
comparison with the existing DNS data resulted in an 
estimate for the constants C BI w 1.46, and c uu ~ 0.36. 



Then, we considered the region ft — > oo. Even though 
such a condition may not be realizable in nature, from 
a methodological point of view, as we will see below, it 
enables to obtain the desired analytical approximation. 
The ft — > oo asymptotic solution with corrections, linear 

in the small parameter ft 4 ^ 3 , were found. Now we are 
armed to suggest an interpolation formula 



Et{ft) 



3/2 



lift 

3 Cut i 



8C R 



^(llft/3c„„) 2/3 +(8C RI )' 

(3.1a) 

that coincides with the exact solutions for ft = and for 
ft — > oo, including the leading corrections to both asymp- 

totics, linear in ft, and ft 4 ^ 3 . Moreover, in the region 
ft ~ 1, Eq. (3.1a) accounts for the structure of the exact 
polynomial. As a result, the interpolation formula (3.1a) 
is close to the numerical solution with deviations smaller 
than 3% in the entire region < ft < oo, see upper mid- 
dle panel on Fig. 1. Together with Eq. (2.17c) it produces 
a solution for £+, that can be written as 



S+(£ + ) ~ (L+)" 



1 + 



(l+/L+f/3) 1 , (3.1b) 



where = 3L+/14 , L\ = 3L+/11 n and k is the von- 
Karman constant. This formula gives the same accuracy 
~ 3%, see upper left panel in Fig. 1. We demonstrate 
below that the proposed interpolation formulae describe 
the ^-dependence of the correlations with a very reason- 
able accuracy, about 10%, for any < ft < oo, see black 
dashed lines in Figs. 1. 

Unfortunately, a direct substitution of the interpola- 
tion formula (3.1) into the exact relation for S e obtained 
from the system (2.17) works well only for small ft, in 
spite of the fact that the interpolation formula is rather 
accurate in the whole region. We need therefore to de- 
rive an independent interpolation formula for S e . Using 
expansions for small ft <C 1 and large ft » 1 we suggest 



q+ c 



Si Q + 6(c uu a)^S+? 



(3.1c) 



(l + al+/L+f /3 
in which 

c+ __ oi/4 r /r> 1 / 4 p+ 

"-> e ,0 — tnu^ue / ^ri «• > 

S+? = -C ue (2C m -(llC ue -3C slJ )/3S+°°L+)/L+, 
S+°° = -14(C SU - 4C ue /3)/3L+ 

and a satisfies 

Stj + = S+° C L+ + SS+^L+icuua)^ 3 - 4aS e , /3 , 

with 

StJ + = -Cue (3/4C M - 22 + 3C SU /C„ e ) /24C RI . 

Equation (3.1c) is constructed such that the leading and 
sub-leading asymptotics for small and large ft coincide 
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with the first two terms in the exact expansions at "al- 
most" neural stratification and extremely strong stratifi- 
cation. As a result, Eq. (3.1c) approximates the exact so- 
lution with errors smaller then 5% for ft < 1 and ft > 50 
and with errors smaller than 10% for any ft, see lower 
left panel in Figure 1. 

Substituting the approximate Eqs. (3.1) into the exact 
relations (2.17), one gets approximate solutions Eg and 
F+ with errors smaller than 10%, see rightmost panels 
in Figure 1. 

B. Mean velocity and temperature profiles 



This approximation is plotted as a dashed line on Fig- 
ure 2 for comparison. One sees that the approxima- 
tion (3.3b) works much better than the traditional one. 
Thus we suggest Eq. (3.3b) for parameterizing meteoro- 
logical observations. 

The temperature profiles in our approach look similar 
to the velocity ones: they have logarithmic asymptotic 
for £ < L and linear behavior for £ > L. Correspond- 
ingly, they can be fitted by a log-linear approximation, 
like (3.3a), or even better, by an improved version of it, 
like Eq. (3.3b). Clearly, the values of constants will be 
different: k k t , Li £i,t, etc. 



In principle, integrating the mean shear Si and the 
mean temperature gradient S£, one can find the mean 
velocity and temperature profiles. Unfortunately, to do 
so we need to know S+ and S£ as functions of the eleva- 
tion z, while in our approach they are found as functions 
of £/L. Remember, that the external parameter £ is the 
outer scale of turbulence that depends on the elevation 
z. The importance of an accounting for the proper phys- 
ically motivated dependancc of £ on z for an example 
of channels and pipes has been recently shown by L'vov 
et al. (2008). For the problem at hands, we can safely 
take £ — z if z <C L, however when z > L the function 
£(z) is not found theoretically although it was discussed 
phenomcnologically with support of observational, exper- 
imental and numerical data. It is traditionally believed 
that for z > L the scale £ saturates at some level of order 
L [see, e.g. Eq. (0.4)]. 

The resulting plots of U + are shown on Figure 2, left 
panel. Even taking £{z) — z one gets a very similar 
velocity profile, see Figure 2, right panel. With £(z) = z 
we found an analytical expression for the mean-velocity 
profile using the interpolation Eq. (3.1b) for 5+: 



U+(z) = - In 



(l + ^/l + (z/L 2 f 3 



- . (3.2) 



Here z u0 is the roughness length. 

The resulting mean velocity profiles have logarithmic 
asymptotic for z < L and a linear behavior for z > L in 
agreement with meteorological observations. Usually the 
observations are parameterized by a so-called log-linear 
approximation (Monin and Obukhov, 1954): 



U + = k 1 ln(z/z u0 ) + z/Li , 



(3.3a) 



which is plotted in Figure 2 by dotted lines. One sees 
some deviation in the region of intermediate z. The rea- 
son is that the real profile [see, e.g. Eq. (3.2)] has a 
logarithmic term that saturates for z>L, while in the 
approximation (3.3a) this term continues to grow. To fix 
this one can use Eq. (3.2) (with L 2 = L\ for simplicity), 
or even its simplified version 



1 



U + = - In 



k Zuo^/l + iz/Ltf 



Li 



(3.3b) 



C. Profiles of second-order correlations 

The computed profiles of the turbulent kinetic and 
temperature energies, horizontal thermal flux profile and 
the anisotropy profiles are shown on Figure 1 in the mid- 
dle and right panels. The anisotropy profiles, lower mid- 
dle panel, saturate at £/L « 2, therefore they are not 
sensitive to the z-dependence of £(z); even quantitatively 
one can think of these profiles as if they were plotted as 
a function of z/L. 

Another issue is the profiles of (upper middle 
panel) and of E e and _F+ (rightmost panels), that are 
oc (£/L) 2 / 3 for £ » L (if realizable). With the interpola- 
tion formula (0.4) the profiles of the second order corre- 
lations have to saturate at levels corresponding to ft = 1 . 
This sensitivity to the z-dependence of £{z) makes a com- 
parison of the prediction with experimental data very de- 
sirable. 



D. Turbulent transport, Richardson and Prandtl 
numbers 



In our notations the turbulent viscosity and thermal 
conductivity, turbulent Prandtl number, the gradient- 
and flux-Richardson numbers are 



7~xz 

v t = 

<~>u 
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(3.4c) 
(3.4d) 
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S l L+Sf si 2 ' 
0F Z J_ _ ft_ 



Rigrad — Riflux P^T 



(3.4f) 



With Eqs. (3.4a) and (3.4b) we introduce also two di- 
mensionless functions C v {ft) and C x {ft) that are taken 
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FIG. 1: Color online. Log-log plots of the normalized velocity mean shears Si = £ + S+ and L+S+ (left upper panel) , normalized 

mean-temperature gradients S| = £ + S& and L + S& (left lower panel), the turbulent kinetic energy E£ and E^/ft 2 3 (middle 

upper panel), partial kinetic energies tu/E k (middle lower panel), temperature energy E+ and E+/ft 2/3 (right upper panel) 

and horizontal thermal flux F£ and F^/ft 2 ^ 3 (right lower panel) vs. ft = £/L — £ + /L + . Red and blue solid lines - exact 
numerical solutions before and after normalization by the large ft asymptotics, black dashed and dot-dashed lines - approximate 
analytical solutions. The region i > L may not be realized in the Nature. In this case it has only methodological character. 




FIG. 2: Computed with Eq. (0.4) (for di = d 2 = 1) plots of U + (blue solid lines) vs \n(z/L) and vs. z/L (inserts) for L + = 1000. 
In the left panel £(z) is taken from Eq. (0.4), while in the right panel £(z) — z. Log-linear approximation (3.3a) to all profiles 
is shown by dotted lines, its improved version (3.3b) by dashed lines. "! =" stands for 7^. The region t^L may not be realized 
in the Nature. In this case it has only methodological character. 
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as ^-independent constants in the down-gradient trans- 
port approximation (0.3) described in the Introduction. 
We will show, however, that these functions have a strong 
dependence on ft, going to zero in the limit ft — > oo as 

l/ft 4 ^ 3 . Therefore this approximation is not valid for 
large ft even qualitatively. 



1. Approximation of down- gradient transport and its 
violation in stably stratified TBL 

As we mentioned in the Introduction, the concept of 
the down-gradient transport assumes that the momen- 
tum and thermal fluxes are proportional to the mean 
velocity and temperature gradients, see Eqs. (0.3): 



-is T S v 



(3.5) 



where v T and Xt are effective turbulent viscosity and 
thermal conductivity, that can be estimated by dimen- 
sional reasoning. Equations (0.3), giving this estimates, 
include additional physical arguments that vertical trans- 
port parameters should be estimated via vertical turbu- 
lent velocity, yjr zz , and characteristic vertical scale of 
turbulence, £ z . The relations between the scales £j in dif- 
ferent j-directions in anisotropic turbulence can be found 
in the approximation of time-isotropy, according to which 



7 



In 



(3.6) 



Here 7 is a characteristic isotropic frequency of turbu- 
lence, that for concreteness can be taken as the kinetic en- 
ergy relaxation frequency 7„ u . The approximation (3.6) 
is supported by experimental data, according to which in 
anisotropic turbulence the ratios £i/£j (i 7^ j) are larger 



then the ratios 



that are close to unity. 



With this approximations v T and % T can be estimated as 
follows: 

V T = C V T ZZ /^ UU , \t = C x T zz /j uu , (3.7) 

where, according to the approximation of down-gradient 
transport, the dimcnsionlcss parameters C v and C' x are 
taken as constants, independent of the level of stratifica- 
tion. 

In order to check how the approximation (3.5), (3.7) 
works in the stratified TBL for both fluxes, one can con- 
sider Eqs. (3.5) as definitions of v T and Xt and Eqs. (3.7) 
as definitions of C v and C v . This gives 



(3.8a) 
(3.8b) 





'Juu 


'fuu 


7~zz 


^7 


T+S+ 


Fx_ 




7+ 

>uu 


7~zz 


So 


TzzSq 



Recall, that in this paper the down-gradient approxima- 
tion is not used at all. Instead, we are using exact bal- 
ance equations for all relevant second order correlations, 



including t xz and F x . Substituting our results in the 
RHS of the definitions (3.8) we can find, how C v and C x 
depend on ft = £/L that determines the level of stratifi- 
cation in our approach. 

The resulting plots of the ratios C„(ft)/C u (0) and 
C x (ft)/ C x (0) are shown in the leftmost panel in Figure 3. 
One sees that the C v {ft) and C x (ft) can be considered 
approximately as constants only for £ < 0.2 L. For larger 
£/L both C v (ft) and C x {ft) rapidly decrease, more or less 
in the same manner, diminishing by an order of magni- 
tude already for £ w 2 L. For larger £/L one can use the 
asymptotic solution according to which 



1 



ft 



2/3 



(3.9) 



£ + , rau — g + — £ + , >zz 

This means that both functions vanish as 

/ T \ 4/3 / T \ 4/3 

a(*i)*0.01^ 7 J , C x (**)~ 0.003 (jj ,(3.10) 

where numerical prefactors account for the accepted val- 
ues of the dimensionless fit parameters. 

The physical reason for the strong dependence of C v 
and C x on stratification is as follows: in the RHS of 
Eq. (2.10e) for the momentum flux and Eq. (2.11b) for 
the vertical heat flux there are two terms. The first ones, 
proportional to t zz and velocity (or temperature) gradi- 
ents correspond to the approximation (3.5), giving (in our 
notations) C v =const and C x =const, in agreement with 
the down-gradient transport concept. However, there are 
second contributions to the vertical momentum flux oc F x 
and to the vertical heat flux, that is proportional to f3Eg . 
In our approach both contributions are negative, giving 
rise to the counter- gradient fluxes. What follows from 
our approach, is that these counter-gradient fluxes can- 
cel (to the leading order) the down-gradient contributions 
in the limit ft — > 00. As a result, in this limit the effec- 
tive turbulent diffusion and thermal conductivity van- 
ish, making the down-gradient approximation for them 
(with constant C v and C x ) irrelevant even qualitatively 
for £ > L. 

In our picture of stable temperature-stratified TBL, 
the turbulence exists at any elevations, where one can 
neglect the Coriolis force. Moreover, the turbulent ki- 
netic and temperature energies increase as (£/L) 2 ^ 3 for 
£ > L, see Figure 1. At the same time, the mean velocity 
and potential temperature change the (£/L)-dependence 
from logarithmic lo linear, see Figure 2 and (modified) 
log-linear interpolation formula (3.3b). Correspondingly, 
the shear of the mean velocity and the mean temperature 
gradient saturate at some elevation (and at some £/L), 
and Rigrad saturates as well. This predictions agree with 
large eddy simulation by Zilitinkevich and Esau (2006), 
where Ri gra d can be considered as saturating around 0.4 
for z/Lk 100. 

Notice that the turbulent closures of kind used above 
cannot be applied for strongly stratified flows with 
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FIG. 3: Color online. Log-log plots of "down-gradient coefficient-functions" C„ (solid blue lines) and C x (red dashed lines) - left 
panels; turbulent Prandtl number Pr T (green lines on middle panels) and Rifl ux (black solid lines), Ri gra d (black dashed lines) 
- on right panels as function of £* — £/L (upper panels) and vs. Ri gra d (lower panels). Notice, that the presented dependencies 
have qualitative character, and the choice of constants C... depends on the actual functional form I (z). For simplicity, we took 
I (z) — z.The region £ > L may not be realized in the Nature. In this case it has only methodological character. 



Rigrad ^ 1 (may be even at Ri gra d ~ 1)- There are 
two reasons for that. The first one was mentioned in 
the Introduction. Namely, for Ri gra d ^ 1 the Brunt- 

Vaisala, frequency N = \/pS , N + = ij S£/L+, is larger 
then the eddy-turnover frequency and therefore there 
are weakly decaying Kelvin-Helmoholtz internal gravity 
waves which, generally speaking, have to be accounted 
for in the momentum and energy balance equations. 



The second reason, that makes the results very sensi- 
tive to the contribution of internal waves follows from the 
fact that vortical turbulent fluxes vanish (at fixed veloc- 
ity and temperature gradients). Therefore even relatively 
small contributions of different nature to the momentum 
and thermal fluxes may be important. 



The final conclusion is that the TBL modeling at large 
level of stratification requires an accounting for turbu- 
lence of the internal waves together with the vortical 
turbulence. Definitely, new observations, laboratory and 
numerical experiments with control of internal wave ac- 
tivity are very likely. 
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APPENDIX A: ON THE CLOSURE PROBLEM 
OF TRIPLE CORRELATIONS VIA SECOND 
ORDER CORRELATIONS 

Let us look more carefully at the approximation (2.12), 
which is 

luu — Cuu V ^-^k I ^ , Tri = Cm'fuu ? (Al) 

7bj = Cri7ri > lee = Ceeiuu , 7rd = C u g^ uu . 

The dimensional reasoning that leads to this approxima- 
tion is questionable for problems having a dimensionless 
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parameter ft . Generally speaking, all "constants" c... and 
C... in Eq. (Al) can be any functions of ft. Presently we 
just hope that a possible ft dependence of these functions 
is relatively weak and does not affect the qualitative pic- 
ture of the phenomenon. 

Moreover, even the assumption (2.5a) that the dissipa- 
tion of the thermal flux ej is proportional to the thermal 
flux and the assumption (2.5b) that the dissipation of 
Eg, e cx Eg are also questionable. Formally speaking, one 
cannot guarantee that the triple cross-correlator (6uu) + 
that estimates e + , can be (roughly speaking) decomposed 
like (u9) sj (uu), i.e really proportional to F = (uff) as it 
stated in Eq. (2.5a). Theoretically, one cannot exclude 
the decomposition (Quu) <~ (uu) (09), i.e. a contri- 
bution to e cx E K . Similarly, the dissipation e in the 
balance (1.5c) of Eg, that is determined by the correla- 



tor (2.5b), is cx (OOu), as it follows from the decompo- 
sition (69u) <~ (99) \J (uu) and is stated in Eq. (2.5b). 
This correlator allow s, for example, the decomposition 
(99u) ~ (9u) y/(96), i.e. contribution to e cx F. This 
discussion demonstrates, that the situation with the dis- 
sipation rates is not so simple, as one may think and thus 
requires careful theoretical analysis that is in our agenda 
for future work. Our preliminary analysis of this problem 
shows that all fitting constants are indeed functions of ft . 
Fortunately, they vary within finite limits in the entire 
interval < ft < oo. Therefore we propose that the ap- 
proximations used in this paper preserve the qualitative 
picture of the phenomenon. Once again, the traditional 
down-gradient approximation does not work even qual- 
itatively because corresponding "constants" C„ and C x 
vanish in the limit ft — > oo. 
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